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Using Brownian dynamics simulation, we study the orientational order in the melting transition of 
. a colloidal system with a 'soft' Yukawa potential. The bond-orientational order parameter $6 and the 

' bond-orientational order function gsir) are calculated in two-dimensional systems. We found that a 

I two-stage transition and a hexatic phase indeed exist in two-dimensional melting, which is consistent 

^N) ■ with the prediction of the Kosterlitz-Thouless-Halperin-Nelson- Young theory. For comparing with 

the melting process in three-dimensional systems, we introduce the probability distribution of the 
O , single-particle local-order parameter. Based on extensive simulations, the breakdown of local order 

• is qualitatively suggested only to occur on the fractional part of the colloidal system for the two- 

dimensional melting, but in three-dimensional melting, that breakdown takes place on the whole 
system at the same time. 
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The melting transition in two-dimensional (2D) systems has received much attention during the past ten years. 
Unlike higher-dimensional systems, 2D solids possess a quasi-long-range positional order. This property implies that 
the melting transition in two dimensions is different from that in 3D systems. In 1973, Kosterlitz and Thouless sug- 
gested that the melting transition in 2D crystals was a continuous process mediated by the dissociation of dislocation 
pairs [1]. Later, the resulting phase of this continuous transition was shown not to be an isotropic phase because it 
G ' still had a quasi-long-range orientational order [2]. In 1979, Young pointed out that a second-transition, which was 
induced by the formation of disclinations, would drive this so-called hexatic phase into liquid [3]. This theory, which 
is well-known as the Kosterlitz-Thouless-Halperin- Nelson- Young (KTHNY) theory [4], thus proposes the following 
O two-stage scenario for melting: The solid first undergoes a continuous transition to become a hexatic phase with a 
quasi-long-range orientational order; then, another continuous transition drives the hexatic phase to a disorder liquid. 
During the years, a large number of experiments and simulations have been performed to verify the KTHNY theory 
PsJ I [5~8]; however, results remain controversial, which It may be caused by the size effect of the systems, because long- 
^ • wavelength fluctuations play an important role in the KTHNY theory, they are cut off in finite-sized simulations. In 
On [ particular, very large-scale simulations of Lennard- Jones systems seem to provide evidence for the existence of the 
. hexatic phase [9]. However it was argued that this transition might depend on the specific properties, such as the 
' interparticle potential of the studied systems [10,11]. Recently, a series of experiments was performed to calculate 
: the elastic moduli and the dislocation interaction of 2D colloidal crystals [12-15]. The renormalized Young's modulus 
\^ ' of the crystals, Kji, was found to be consistent with the KTHNY theory, and the dissociation of dislocations was 
observed experimentally. 

' This work is devoted to a study of the melting transitions of charged colloidal crystals. Colloidal crystals are known 
] to be ideal model systems for statistical physics. A particular advantage of such systems is that, due to the particle 
H ■ size, the individual colloid positions and motions can be directly observed. Moreover, the inter-particle potentials 
in such colloidal crystals are in most cases precisely known. For example, the well-known Dejaguin-Landu-Verwey- 
Overbeek (DLVO) potential gives a good description for the effective pair interaction of one-component model of the 
colloidal systems. However, for a dilute charged-stabilized colloidal system in which many-body interactions can be 
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Q , ignored, a pairwise Yukawa potential, which is only the electrostatic part of the DLVO potential, is suitable [16-19]: 

V{r) = Uoaeyip[~X{r-a)/a]/r, (1) 



X 



5_j • where Uq sets the energy and a the length scale. The screening parameter A describes the 'softening' of the particles: 
Ci When A increases from zero to infinity, the cores of interacting particles change from very soft to extremely hard [20] . 

While most of the previous research focused on hard-core colloidal particles, we study a 'soft' disk system with 
Yukawa interactions. However, no consensus was achieved even in such systems [10,21]. In this paper, the bond- 
orientational order parameter $6; the probability distribution of the single-particle order parameter P(0e), and the 
orientational correlation function gsir) are calculated to describe the melting processes of 2D 'soft' colloidal crystals. 
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The translational order is ignored because it cannot describe the properties of a two-stage melting transition. Melting 
in 3D systems was also studied for comparison with 2D systems. 

First, we briefly describe the standard Brownian-dynamics simulation method [20,22,23]. It is based on a finite- 
difference integration of the irreversible Langevin equations. If the effect of the hydrodynamic interactions is ignored, 
the equations of motion of N particles at equilibrium are derived from Eq. (1) as 

^Vi{t) = F,{t)+-R{t), (2) 

where i = 1, ....iV labels the A'' particles, ^ is the friction coefficient, R(t) is the Langvein random force of the solvent, 
and Fj is the total inter-particle force on particle i. The finite diff'erence integration of Eq. (2) is 

U{t + At) = v,{t) + (l/e)F,(f)At + (Ar)fi + 0{{Atf), (3) 

where (Ar)/{ is a random displacement sampled from a Gaussian distribution of zero mean and variance 



(Ar)| = 4Z?oAi. (4) 

Do = ksT/^ is the short-time diffusion coefficient, where /cg is the Boltzmann constant and T is the temperature. 

In all simulations, we used reduced units such that a = 1,Uq = 1, and p = N/V = 1. The screening parameter A = 8. 
In the case of 2D simulation, we used a periodically repeated rectangular simulation box of volum e V with N = 2500 
particles and started from a triangular lattice; The cutoff Vc was set as 4. If, where f = yJV/N = p~'^l'^ . The scale 
of the simulation box was set to x : y = \f2> : 2 in order to minimize the efi'ect of the boundary. Particle trajectories 
were generated according to Eq. (3), and the simulation data were gathered within the range of 5tb (jb = <j^^/Uq) 
after a long equilibration time (more than 45tb). In our simulations, we varied the reduced temperature T* while the 
other parameters a, Uq, p, and A were fixed. Several different runs corresponding to different reduced temperatures 
were done. 

There are some remarks to be made on the simulation details. First, the timestep At depends on the screening 
parameter A. When the potential becomes 'hard', a shorter At must be chosen. Another important thing is the 
relaxation time of the 2D system. It should be very long so that equilibrium will be reached if a two-stage continuous 
melting transition occurs. In particular, the 'hysteresis loops' will be incorrect when the system is heated or cooled 
too fast [24]. Lastly, the center of mass of the system was fixed during the simulation in order to avoid super drift. 

Melting in 3D systems was also simulated for comparison with 2D melting. The method is the same as in 2D 
systems except that the 4 on the right-hand side of Eq. (4) is changed to 6 and f = p~^^^ in the cutoff length 
Vc = 4. If. We used a system of 500 particles with periodic boundary conditions and started from FCC structures. 

The bond-orientational function gsir) is used to identify the existence of the hexatic phase [25]. It is defined as 



(V6*(r')V6(r'-r)) 
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where g{r) — {S{r')5{r' — r)) is the pair distribution function and tpei^) is the local bond orientational order parameter 
and is given by 



Here, Nh denotes the number of nearest neighbors of the mth particle, and 6mn is the angle between some fixed axis 
and the bond joining the mth particle with an nth particle. According to the KTHNY theory, the bond-orientational 
function gsir) will have an algebraic decay in a hexatic phase and an exponential decay in the liquid phase. 

The order parameter $6 was introduced by Nelson and Halperin to characterize the structural order in 2D systems 
[2]. It is given by [26] 

\ m=l ° n=l / 

When the system belongs to the fiuid phase, |$6|^ ^ 1- On the other hand, |$6|^ ^ 1 means that the system is a 
perfect crystal. Especially, we define a single-particle order parameter (pe, which is the ensemble average of the local 
bond-orientational parameter ipQi 
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Clearly, 0g is the order parameter of a single particle, and one can calculate 4>e for each particle in our simulation 
box. 

As $6 in two dimensions, the bond order parameter Qq is used to describe the orientational order of a 3D system. 
It is defined by calculating the spherical harmonics over all bonds [27]. Value of Qe from 0.354 to 0.663 represent 
systems with crystal structures [28]. For example, when the system is a perfect FCC crystal, = 0.57452, and when 
system is a liquid, Qe <^ I. As 06 in 2D systems, we also introduce a single-particle local order parameter qq, which 
is the order parameter of a single particle. The probability distributions of gg and 06 represent a possible way to 
qualitatively study 2D and 3D melting transitions, as shown in the following context. 




FIG. 1. Bond-orientational function ps(r) versus r/a for various values of the reduced temperature T* = ksT /Uo. The 2D 
simulation system consists of 2500 particles and the screening parameter A = 8. The straight line, with slope —1/4, is a guide 
for the eyes. 



The 2D simulation results for the bond-orientational functions gsif) are given in Fig. 1 with different reduced 
temperatures T* — kBT/Uo. At the low reduced temperature, the bond-orientational order function gsir) does not 
decay, which implies the system is in a solid phase with a long-range bond-orientational order, but when the reduced 
temperature rises to 0.600, gsii") decays algebraically with an exponent less than —1/4, the upper limit of the existence 
of a stable hexatic phase predicted by the KTNHY theory. This provides very strong evidence for the hexatic phase. 
As the temperature keeps increasing, the system becomes a disordered liquid, and (7_b(?') decays exponentially. The 
exponent of the algebraic decay always decreases with increasing temperatures in the hexatic phase. However, this 
decrease of the exponent is very small because the temperature region of the haexatic phase is only 0.008. As a result, 
the temperature dependence is hardly observed in the simulation results. It should be noted, in Fig. 1, that gsir) 
becomes flat for larger r due to the finite-size effect, such as gsir) at the hexatic phase for r > 20.00. In Fig. 1, we 
only present gsir) for r < 22.00. The length of the x-axis of our simulation box is 55.836 x a. 

Note that because of the large fluctuation near the critical point, the phase boundaries we obtained are not very 
precise, but what we are interested is the process of the melting transition. Additionally, we perform extensive 
simulations in very large ranges of T* from 0.050 to 1.200 and of p from 0.6 to 1.2. Moreover, a two-stage melting 
transition is clearly observed in each case. Thus, we can conclude that melting transition of a soft disk system is 
absolutely a solid-hexatic-liquid process. 

It is worthwhile to note that the gsir) decays algebraically with an exponent is —2/3 at T* — 0.609, which is faster 
than the —1/4, predicted by the KTHNY theory. In our simulations at this reduced temperature, it took a very long 
time for the systems to reach equilibrium. We conjecture that T* ~ 0.609 is the critical reduced temperature for the 
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transition from the hexatic phase to a disordered phase. The detailed properties of the system at T* = 0.609 or the 
hexatic-disordered transition are still unclear and need to be studied further. 




FIG. 2. (a) Probability distributions of the single-particle order parameter 4>(i in 2D melting. P{4>(i) exhibit an orderless 
distribution at T* = 0.608, which implies the existence of a hexatic phase at this temperature. The inset shows the value of the 
order parameter |$6p at different T* . Two drops occurs, one at T* — 0.600 and the other at T* = 0.608, implying a two-stage 
transition, (b) Probability distributions of qe in 3D systems. The values are distributed over a very narrow range. The inset 
shows the value of Qe at different T*. The collapse at T* — 0.440 represents the occurrence of a first-order transition. 



In Fig. 2, we present the simulation results for the bond-orientational order parameters in melting of 2D and 3D 
colloidal crystals. At low temperature, |<i>6p has a large value, which implies that the 2D system is in a solid phase [see 
the inset in Fig. 2 (a)], and it decreases slowly with increasing T*. The first drop is found at T* = 0.600, corresponding 
to the fact that a transition occurs at this temperature. Then, |$6p returns to a slow decrease until the second drop 
takes place at T* — 0.608. We suggest that a two-stage phase transition, indeed, occurs in the melting of 2D systems 
with soft Yukawa potentials. In contrast, the single collapse of Qe at T* — 0.440, as shown in the inset of Fig. 2(b), 
implies that 3D melting with a soft Yukawa interaction is just a well-understood first-order transition. Obviously, the 
bond-orientational order parameter is very sensitive to the structure transition of the simulation systems, which can 
be used to confirm the emergence of two-stage melting. 

The probability distribution of the single-particle order parameter P{(I)q), as shown in Fig. 2(a), exhibits a Gaussian- 
type behavior of both high and low temperatures, and the variance increases with temperature. However, when the 
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system is at an intermediate T*, P{cf)(i) is not a Gaussian- like curve anymore. We conjecture that the orderless distri- 
bution in Fig. 2(a) is caused by the appearance of isolated dischnations (5— and 7— fold coordinated particles having 
exclusively 6- fold coordinated neighbors), which is characteristic of the hexatic phase destroying the translational 
symmetry of the crystal, as observed in Ref. [1]3. In the case of 3D systems, P{q6) keeps a Gaussian-like distribution 
for the entire melting process [see Fig. 2(b)]. This appears to be due to an enhancement of the thermal motions of 
all the particles with increasing temperatures. Also the breakdown of the ordered phase takes place for the whole 
systems at the same time. Moreover, it is worth mentioning that, in all the simulation results, the variance of the 
distribution of melting in a 3D system is obviously smaller than that in a 2D system. Additionally, we suggest that 
P{4>6) is more intuitive and easier than pB(r) for presenting detailed structural information in the melting process. 

In conclusion, we have performed Brownian dynamics simulations on the melting of 2D colloidal crystals in which 
particles interact with a 'soft' Yukawa potential. A two-stage melting transition and a metastable hexatic phase 
between the solid and the liquid phases arc observed in our simulations, which are consistent with the KTHNY 
theory. On the other hand, the melting in 3D soft Yukawa systems is just a one-step phase transition. In addition, 
through simulations of the bond-orientational order parameter and the probability distribution of the single-particle 
local-order parameter, the 3D melting is thought to occur with the breakdown of the locally ordered phase over the 
entire system at the same time. 
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